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Abstract 

Full-potential LMTO calculations for KNVJO3 reported up to now provided 
a reasonable description of off-center equilibrium displacements, including 
transversal optical T phonons. However, when addressing more sensitive 
phonon properties and the behavior of the ferroelectric instability over the 
Brillouin zone, the need was realized to achieve a substantially higher level of 
accuracy. 

In order to arrive at ultimately accurate results available with the LMTO 
method in the local density approximation, the stability of full-potential 
LMTO predictions for off-center displacements in KNb03, as depending on 
the choice of basis and expansion cutoffs, has been thoroughly investigated. 

With the calculation setup thus optimized, supercell frozen phonon calcu- 
lations aimed at the study of the chain-structure instability over the Brillouin 
zone have been done, and the long-wavelength limit of the LO phonon is 
discussed. 

I. INTRODUCTION 

In referring to the nature of the ferroelectric instability in KNbOs, which develops as 
Nb atoms are displaced from their symmetric positions in the centre of the cubic perovskite 
cell, it became usual to mention a crossover from displacive to order-disorder type of tran- 
sition, that gave rise to considerable controversy ^0,0,^,0]. Indeed, both the temperature 
dependence of the dielectric constant and the anisotropy of the transverse acoustic and 
lowest optical phonon branches may be understood in the most natural way in the frame- 
work of the displacive model 0. At the same time, recent stimulated Raman scattering 
experiments [|§|P|,|TU| clearly suggest that separate off-center potential minima for Nb atoms 
exist, and the hoppings between them can be induced by an appropriate photon pump- 
ing. This controversy may be understood in the following way: The displacive scenario 
supposes that the (in our case Nb) atoms are sitting at, or rather oscillating about, the 
central symmetric position in the high-temperature paraelectric phase, and this is definitely 
not the case in KNbOs. In fact Nb atoms are displaced in each individual unit cell all 
way around through the sequence of phase transitions. But the textbook description of 
the order- disorder model with Nb atoms statistically distributed over eight [lll]-displaced 
positions is also not true in its implicit assumption that each individual atom has a choice 
of eight local potential minima around the centre of the perovskite cell. Be it so, the scan 
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of the potential hypersurface related to the displacement of a single individual atom in an 
otherwise non-polar paraelectric crystal would clearly exhibit these minima. Such scanning, 
that is difficult to implement experimentally but easy to imitate in a supercell total-energy 
calculation, shows no traces of off-center minima and no noticeable anharmonic damping 
of local potential surface related to a single atom displacement beyond any uncertainty of 
calculation ||11|| . At the same time, precise total-energy calculations for perfectly ordered 
ferroelectric KNbC>3 [I2], 13| , |n| routinely reproduce the tendency of Nb atoms to displace 



simultaneously off-center in the [lll]-direction, and the T-TO frequencies which come out 
reasonably well in such calculations [|T^,|T^,|T^,|r^] suggest that the potential hypersurface 
related to various atomic displacements is described sufficiently well. 

The message from this mixed evidence is that the off-center potential minimum for 
the Nb atom in each individual perovskite cell is determined by the electric field in some 
vicinity of the cell in question. In other words, the displacement of Nb atoms is correlated, 
and the ferroelectric instability develops as the softening of particular phonon modes. This 
statement belongs to the vocabulary of the displacive model. The important difference is that 
the particular vibration pattern may have finite spatial extent, or correlation length, even at 
very low temperatures. Another difference is that a particular soft mode may not necessarily 
freeze down with temperature, but rather remain frozen at all temperatures throughout the 
ferroelectric transitions. What evolves with temperature is either the statistical mixture, 
or the dynamical interplay, of different frozen modes, apparently accompanied by a gradual 
decrease of the correlation length. 

With respect to spatially correlated displaced Nb atoms, a guess that they may build 
chain-like structures has been put forward long ago by Comes et al. |18|] and since then 
addressed by new measurements (see, e.g., Ref. [|19|)- It was not however possible to extract 
experimentally any information as for the correlation length along the chain, nor to whether 
there is any correlation between different chains. From the lattice-dynamical point of view, 
the study of the spatial correlation along the chain has to do with the instability region, i.e. 
that of the soft-mode freezing, in some vicinity of the Brillouin zone center. Yu and Krakauer 
|20| mapped the phonon spectrum of cubic KNb03 over the whole Brillouin zone in an ab 
initio calculation using a linear response approach. Once the general shape of the instability 
region is known, additional information on the correlation in the atomic displacements, 
which may as well include the effects beyond the harmonic approximation, may be obtained 
from supercell calculations. In the present paper, we provide a preliminary attempt of such 
analysis, based on the use of a fast full-potential linear muffin-tin orbital (LMTO) method. 



II. REFINED CALCULATION SETUP FOR THE FULL-POTENTIAL LMTO 

CALCULATIONS OF KNbQ 3 



Using the fast and accurate full-potential LMTO code by M.Methfessel [pT] , f22"f as a tool 
for ab initio total energy calculations, we were able to study earlier in a series of papers 
lT].p~3|.p~5|,[T6[1 different aspects of ferroelectricity in KNbOs, including magnitudes of equi- 



librium off-center atomic displacements, coupling of the displacements to the (tetragonal) 
lattice strain, and T-TO phonon frequencies in different phases of this compound. Keeping 
in mind the planned applications to large supercells, we tried to analyze the resources of the 
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method and of our previously used setup concerning an eventual decrease in the computa- 
tional demand without substantial loss of accuracy. At the same time, we tried to find a way 
to improve the description of the 0-K bonding, which came out somehow too rigid in our 
previous calculations. As a result, the coupling of K to the oxygen cage on a ferroelectric 



transition was overestimated - see the discussion in Ref. Ju^ . 

We considered the refinement of our calculation setup along the following guidelines: The 
choice of the muffin-tin spheres sizes is based on the spatial distribution of the potential over 
the cell (with the exception of K whose sphere is enlarged in order to include some interstitial 
space around it) and need not to be changed compared to as discussed in Ref. |TJ|; the same 
applies to the kinetic energies of the Hankel function envelopes (—0.01, —1.0 and —2.3 Ry) 
which provide sufficient variational freedom. As has been discussed at length in Section III of 
Ref. [|13| , one possible source of uncertainty in the calculation results is the termination of the 
spherical harmonic expansions of the charge density in the interstitial region at some value 
L max (in principle different for each inequivalent atom type). The usual way to control this 
is to check the convergency of results depending on this parameter. Another option in the 
adjustment of the calculation setup is the choice of the basis LMTO's. As a reference point 
representing convergency with respect to both these factors, we considered the calculation 
with L max =6 on all atoms, and all valence-state basis functions (s, p and d) included in the 
basis set for each value of the Hankel function energy. The systematic analysis of the trends 
in the calculated properties (we concentrated mostly on the magnitude of the equilibrium 
off-center Nb displacement from the cubic phase and the corresponding energy lowering) has 
then been done in order to find a more economic setup, with respect to L max and the basis 
choice, that does not result in a noticeably loss of accuracy. 

For the K atom with its uncommonly large muffin-tin sphere, we found it preferable to 
keep L max =6, whereas for Nb and O atoms with much smaller spheres the cutoff L max =4 
suffices, as is consistent with the experience of other calculations performed with the code 
by M.Methfessel. As regards the basis set, we found it necessary to include all 5s, 4p and 
Ad states combined with all three Hankel function tails on the Nb site, i.e. adding p and d 
states matched with the Hankel function of the lowest tail energy (—2.3 Ry) as compared 



with the prescription of Ref. [13|]. The choice of basis functions related to K and O sites 



did not need to be changed. Note that the inclusion of the 3d states on O sites for at least 
one (the highest) envelope function energy is essential for describing the very existence of 
the ferroelectric instability. We checked that the further enhancement of the O basis set 
using the 3d functions matched to other Hankel envelopes does not lead to any substantial 
changes. On the K site, the basis functions matched with the lowest-lying Hankel envelope 
may be thrown out whatsoever. 

These refinements do not affect the total energy dependence on volume in any noticeable 
way, nor do they result in substantial changes of the total energy vs. Nb displacement 
behavior in cubic or tetragonal phases of KNb0 3 , compared to the results of Refs. |]i3| . |i"5| . |TTf . 
However, for the orthorhombic phase characterized by much smaller energy differences on 
displacement and tiny magnitudes of the equilibrium displacements, our present results may 
be considered as some improvement over those preliminarily reported in Ref. [|l(| . 
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TABLE I. Positions of atoms in orthorhombic phase of KNDO3 (in terms of lattice parameters) 
as determined by neutron diffraction measurements, Ref. [23|, and optimized in the FP-LMTO 
total-energy calculation. 



Atom 



A 



exp 



A, 



calc 



1 

2 
3 
4 
5 



K 

Nb 

0(1) 

0(11) 

0(11) 





1 

2 



1 

2 

1 

2 







l + A y 



1 
2 

| + A, 
| + A 2 



A, 
A, 



0.0138 

0.0364 
0.0342 
-0.0024 



0.0314 

0.0351 
0.0353 
-0.0037 



III. EQUILIBRIUM STRUCTURE OF THE ORTHORHOBMIC PHASE OF 

KNb0 3 AND T-TO PHONONS 

With the calculation setup as described above, we performed the total energy-driven 
optimization of atomic positions in the orthorhombic phase. In the setting used by Hewat 
p3| , the three lattice vectors are a=3.973A along x=[100] of the cubic aristotype; 6=5.695A 



along y=[011], and c=5.721A along z=[011]. The positions of the atoms have been optimized 
in a four- dimensional parameter space as done earlier in Ref. i.e. keeping the lattice 
volume and strain constant in agreement with the experimental data. The displacements 
of atoms from their symmetry positions in terms of corresponding lattice vectors are shown 
in Table |. As compared with earlier FP-LMTO calculations with a slightly different setup 
]Tj| , the ^-displacement of K was somehow reduced and the z-displacement of both 0(1) 
and O(II) increased, resulting in better agreement with the experimental geometry; the 
estimate of the small O(II) displacement in the y direction also became more accurate. Of 
all structure parameters, only our estimate of the K ^-displacement still remains beyond the 



experimental uncertainty as given in Ref. f23 



Based on this refined geometry, we calculated the T-TO phonon frequencies as another 
benchmark for the quality of our description of the total energy hypersurface as a function of 
atomic displacements. The symmetry coordinates for the displacements from the equilibrium 
positions in the Amm2 space group of orthorhombic KNbOa are as given in Ref. [|16[]. Based 
on a polynomial fit of the total energy hypersurface in terms of symmetry coordinates, we 
calculated the frequencies of the phonon modes compatible with the B 2 (that including the 
soft mode) irreducible representation within our new calculation setup. The frequencies and 
eigenvectors are shown in Table O. As compared to our earlier calculation ||16|| , we obtained 
somehow better agreement with the experimental frequencies. Moreover, the contribution of 
the K displacement in the eigenvector of the soft mode decreases, which means that the soft 
mode represents primarily the vibration of Nb against the oxygen sublattice, with relatively 
unaffected K in the background. This seems to be in better agreement with the displacement 
pattern of atoms in the rhombohedral structure that emerges as the B 2 soft mode freezes 
down, and thus fixes the unaccuracy pointed out in our previous study where the coupling 
of K to the oxygen octahedra was overestimated (see Fig.l of Ref. [(Uj). The calculated 
frequency of the A 2 mode is 252 cm -1 , which is again in somehow better agreement with 
the experimental estimate of 283 cm^ 1 |24| than 224 cm -1 as calculated in Ref. |16 . 
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TABLE II. Calculated T-TO frequencies and eigenvectors for the B<i modes in orthorhombic 
KNb0 3 . 





Eigenvectors (present 


work) 




u calc. (cm 1 ) 


uj exp. (cm 1 ) 


K 


Nb 


Oi 


o 2 


o 3 


0.14 


-0.62 


0.63 


0.32 


0.32 


169* 


soft 


-0.82 


0.32 


0.47 


0.02 


0.02 


204 


195 a -197 c 


-0.30 


-0.05 


-0.54 


0.56 


0.56 


607 


511 a -516 c 



Infrared spectroscopy; Ref. 



IV. CHAIN STRUCTURES AND OFF-r FROZEN PHONONS 

A complete mapping of the phonon frequencies over the Brillouin zone of cubic KNb0 3 
has been done by Yu and Krakauer J20| in an ab initio linear response calculation. Straight- 



forward supercell calculations are too time consuming for any dense mapping of this kind, 
nor are they feasible for an arbitrary k point. They may be however useful for going beyond 
the harmonic approximation for phonons, which is a natural limitation of a linear response 
approach, and for analyzing the ferroelectric instability in more detail for several selected 
k points. In principle, the continous mapping over the Brillouin zone is as well possible 
based on a selection of representative supercell calculations for the fitting of the total energy 



expansion coefficients, as has been proposed by Wei and Chou |25| . This is only a question 
of the computational effort, but not a principal limitation. 

In the present paper, we provide some preliminary results concerning the chain-structure 
instability in orthorhombic KNb0 3 . We performed the calculations for the experimental 
lattice parameters as discussed in Section |TJ, i.e. for the cell volume very close to that used 



in Ref. |[20| . Since the orthorhombic distortion is rather small, one should not expect large 
differences in what regards the region of ferroelectric instability. We concentrated therefore 
on the total energy lowering as a function of finite atomic displacements, with the ultimate 
aim to study the effect of displacement correlations on the anharmonic stabilization of 
structure. For the construction of supercells and the underlying symmetry analysis, we found 
the software tools developed by Stokes |2B| very useful. When considering the frozen phonon 
for a k point along a [100] or [010] symmetry line in an orthorhombic perovskite structure, the 
soft mode is mixed among 9 coupled symmetry coordinates, that demands a corresponding 
fitting of the total energy hypersurface. For preliminary estimates, we considered only [100]- 
displacement of Nb, that dominate in the soft-mode eigenvector as is known from the single- 
cell calculations. The displacements of Nb atoms within particular patterns, depending on 
the given k vector, are from their equilibrium positions in the orthorhombic structure as 
discussed in Section |TT|. 

In order to check how the ferroelectric instability varies over the (100) plane in the 
reciprocal space, we considered the points Y = [010] and S = [0~], that come into M and 
X of the cubic structure in the limit of small orthorhombic distortions. In the direct space, 
corresponding arrangements are infinite [100] chains of stockpiled primitive cells, with Nb 
atoms displaced along [100] in half of the chains and in the opposite direction in the other 
half. With respect to each other, these "up" and "down" chains either pack in the chessboard 
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FIG. 1. Relation between cubic and orthorhombic settings in direct and reciprocal space (top), 
and schematic orientation of dipole moments within chains for supercells corresponding to T, S 
and Y k- vectors (bottom). Note that the electric field along chains is on top of the macroscopic 
field along the z axis. 



configuration (Y), or stick to form interchanging planes (S). The relations between cubic 
and orthorhombic settings and the relevant real-space superstructures are shown in Fig.|I| 
One can expect from the analysis of Ref. 
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that the force constant related to the soft 
mode is not changed significantly along the T—X—M directions. It follows from our earlier 
study of different Nb displacement patterns |16| that the force constants corresponding to 
r and Y are indeed very close 



16 ). However, the 



(see Fig. 3, patterns a and b, of Ref. 
chessboard confguration of chains turns out to be more stable, allowing larger magnitude of 
Nb displacements within the chains. The present study supports this conclusion, now with 
the magnitudes of displacements and the energy gains somehow adjusted with our refined 
calculation setup, and the data for the S frozen phonon added. The total energy gain as a 
function of Nb displacements is shown in Fig. As may be intuitively expected, the curve 
corresponding to the S phonon is the intermediate one between those for T and Y, but 
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Nb displacements (a.u.) 



FIG. 2. Total energy difference (per formula unit) as a function of the [100] Nb displacement 
for the TO T, S and Y frozen phonons. 



rather close to the latter. This type of mutual arrangement can be easily understood from 
electrostatic considerations. Electric dipoles, created (on top of the [001] macroscopic field) 
by additional [100] or [100] Nb displacements in individual orthorhombic cells, tend to form 
chains, but dipoles from different chains, when the latter come closer, have an energy gain if 
oriented antiparallel. The chessboard configuration realizes the maximal separation between 
parallel dipoles on a (100) two-dimensional lattice. Due to the anisotropy of the dipole field, 
the dipole-dipole interaction between chains is much weaker than inside a chain. As a 
consequence, the ideal chessboard arrangement is not bound to occur at temperatures of the 
order 300 K ~1.9 mRy; instead, the chains may be considered as relatively independent (and 
irregular) in their orientation, as long as the number of "up" and "down" chains remains 
roughly equal [^7j . 

This consideration effectively sets the correlation length within the (100) plane, i.e. in 
the direction normal to chains, to one lattice parameter. The analysis of the correlation 
length along chains demands to go away from the (100) plane in the reciprocal space, e.g., 
along M— >R in the cubic setting, and considers the hardening condition of the corresponding 
transversal phonon mode. Based on a uo 2 = condition for a harmonic frequency at a critical 
correlation length, Yu and Krakauer [f2(]] set the latter at about 5 lattice spacings. For the 
orthorhombic lattice, we did the corresponding analysis with a frozen phonon calculation 
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for several q points on a Y— >T line. The values of the force constant related to the [100] 
displacement of Nb as calculated for several increasingly larger supercells of this type are 
shown in Fig.|3] and connected (as a guideline) by a dashed curve labelled Y + k x . The 
supercells have the arrangement of chains as is shown in Fig.|I| in the (100) plane; the [100] 
Nb-displacements within the chains are however now varied as A x (Nb)= ucos(k x x). Each 
supercell represents therefore a body-centered orthorhombic structure with 2, 4 or 8 formula 
units per primitive cell; the force constant is the total energy per unit cell twice differentiated 
in the displacement amplitude u. For k x = 0, the value —65.6 mRy/a.u. 2 obtained from 
the double-cell calculation of the Y phonon (as shown in Fig.|2|) is plotted. The sketches 
in Fig.^| below show the [100] Nb displacement pattern in the supercells. One can see that 
the ferroelectric instability develops at k x < 0.11, i.e. the critical half-period of the Nb 
displacement wave that creates a self-supporting instability is between 4 and 5 unit cells. 
This is in agreement with the estimate done by Yu and Krakauer |[20|| . For a more accurate 
analysis however, one should construct and diagonalize the 9x9 dynamicial matrix and 
search for the condition of a zero eigenvalue to occur. 

As another set of examples for off-r frozen phonons, we considered several A points on a 
[100] line, that corresponds to the uniform Nb displacement pattern over the (100) plane, but 
with harmonic variation (represented by increasingly larger supercells) of A,j.(Nb)~ cos(k x x) 
in the normal direction. The values of the force constant corresponding to this setting are 
shown in FigFJ as laying at a solid line labelled T + k x . Since in this case the phonon 
mode in question is longitudinal, it should not be expected to get soft even at small values 
of k x - The problem of interest in this case is the extrapolation of the force constants 
towards those for the T-LO phonon. The conventional way of treating the T-LO phonons is 
via correcting the dynamical matrix to Born effective charges (see, e.g., [17,28]). Although 



physically transparent, this method incorporates parameters which have to be extracted from 
different schemes. Born effective charges are usually fitted from calculations of the finite 
bulk polarization; the optical macroscopic dielectric function £ oc (q, a;) is either extrapolated 
to oj = from experimental data (Zhong et al. cite the value £oo=4.69 for KNb03 WWli 
or obtained from a linear response calculation, with a typically overestimated, due to the 
local density approximation, value (Yu and Krakauer report £00=6.34 0]). The supercell 
calculation for the LO phonons, on the contrary, involves only the total energy of the ground 
state (in principle, without any limitation as regards the anharmonicity) , and the effects 
of long-range polarization and dielectric response are already implicitly included. Taking 
into account that the force constants should be symmetric with respect to k x — > —k x , we 
extrapolated the force constant related to the [100] Nb displacement in the T-LO mode to 
be 177.2 mRy/a.u. 2 (see Fig.|3]). The corresponding force constant in the T-TO mode (as 
extracted from a single-cell calculation; see Fig.|2]) is —52.3 mRy/a.u. 2 . Based on the relation 



between LO and TO force constants (see Ref. ||29[| , or the relevant formulas (2) and (3) in 
Ref. []3), we arrive at the following relation between the xx-component of the Nb dynamic 
effective charge Z* Nb and e^: (Z* Nb ) 2 /e^ = 3.94. For the two above cited choices of e^, this 
leads to an estimate of the relevant Born effective charge between 4.3 (£00=4.69) and 5.0 
(£oo=6.34). One should note that the Nb-related effective charge tensor is no more diagonal 
in the orthorhombic structure, and we make no estimates for other elements of it than Z xx 
at the moment. Generally, the value of the xx element in the orthorhombic phase may 



be expected to be smaller than 9.13 PS| to 9.37 pi as calculated for the cubic perovskite 
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FIG. 3. Force constant of the [100] Nb displacement for TO [0 1 k x ] (dots on a dashed line) and 
LO [OOfej;] (dots on a solid line) phonons as determined from supercell calculations. The supercells 
used and corresponding Nb displacement patterns are shown below. 
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structure by different methods, because the off-center displacement of Nb from the center 
of a cubic perovskite cell to its position in the orthorhombic cell (along [001], in the present 
setting) was already accompanied by some charge transfer that reduces the polarizability of 
the Nb-0 bonds in perpendicular directions. Moreover, the [100] displacement of Nb in the 
orthorhombic phase is not exactly toward an oxygen atom, as it was the case in the cubic 
structure. 

V. CONCLUSION 

Based on our previous experience in total-energy calculations with the full-potential 
LMTO code by M. Methfessel and on the systematic study of the effects of basis choice 
and the spherical harmonic cutoffs, we attempted to achieve the highest degree of accuracy, 
attainable with the present status of the code, in the description of microscopic structure of 
ferroelectric KNbC>3 in the room-temperature orthorhombic phase. 

This resulted in somehow more accurate, as compared to our previous data, description 
of the equilibrium geometry, i.e. of optimized off-center displacements of atoms, with the 
lattice parameters and the orthorhombic strain kept fixed. The T-TO phonon frequencies 
calculated for the B 2 mode are on average in the same degree of agreement with experimental 
data as in earlier calculations, and the frequency for the A 2 mode has improved. Making use 
of the computational efficiency of the code, we attempted to model several more complicated 
frozen phonon patterns in the supercell calculations with up to 40 atoms, aiming at the study 
of the chain instability and of long- wavelength LO modes. The correlation length for the 
onset of the ferroelectric instability is estimated to be between 4 and 5 lattice constants in 
the [100] direction, i.e. along the chains of displaced Nb atoms, which is consistent with 
the results of a linear response calculation by Yu and Krakauer. The correlations between 
different chains (in the perpendicular direction) are found to be unimportant for the onset 
of the ferroelectric instability, but they stabilize the direction of macroscopic polarization 
compatible with the orthorhombic symmetry, as long as the orthorhombic lattice strain is 
fixed. 

The comparison of force constants, as obtained directly for the T-TO mode and extrap- 
olated to k — > from supercell calculations for a LO mode, makes it possible to exploit 
the relation between Born effective charges and the optical dielectric constant as a tool to 
estimate whichever of these properties is unknown for the system in question. 
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